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Abstract 

Generally  chemical  calculations  are  made  simpler  by  invoking  the  Born- 
Oppenheimer  Approximation,  in  which  the  dynamics  of  electrons  and  nuclei  are 
considered  separable;  this  class  of  chemistry  is  known  as  adiabatic  chemistry.  However, 
in  some  situations  this  approximation  fails  to  effectively  describe  a  chemical  system;  this 
class  of  chemistry  is  known  as  non-adiabatic  chemistry.  Examples  of  non-adiabatic 
chemistry  include  open-shell  reactions  with  atomic  oxygen,  O+N2,  such  as  might  happen 
in  the  upper  atmosphere.  The  B+FT  system,  the  focus  of  this  thesis,  is  also  one  for  which 
non-adiabatic  effects  are  important,  and  was  initially  studied  for  its  possible  use  as  a  High 
Energy  Density  Material  (HEDM). 

The  Hamiltonian  operator  that  describes  chemical  systems  can  be  split  into  the 
sum  of  kinetic  and  potential  energy  operators.  In  order  for  the  Hamiltonian  operator  to  be 
useful  for  creating  solvable  differential  equations  for  the  dynamics  of  a  system,  the 
kinetic  energy  operator  must  be  diagonal.  In  the  adiabatic  representation,  the  potential 
energy  operator  is  diagonal,  but  the  kinetic  energy  operator  is  not.  Chemistry  in  this 
representation  is  only  useful  when  application  of  the  Born-Oppenheimer  Approximation 
allows  the  assumption  that  the  off-diagonal  terms  of  the  kinetic  energy  operator  are 
negligible.  This  assumption  fails  when  the  off-diagonal  terms  of  the  kinetic  energy 
operator,  known  as  non-adiabatic  derivative  coupling  terms  (NAD  terms)  become 
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significant  and  cannot  be  neglected.  This  occurs  when  the  potential  energy  surfaces  of  a 
system  come  close,  touch,  or  even  cross.  In  order  to  form  useable  dynamic  equations 
with  the  Hamiltonian  under  these  circumstances,  it  must  be  represented  in  a  new  basis  in 
which  the  kinetic  energy  operator  is  diagonalized.  Diagonalization  of  the  kinetic  energy 
operator  causes  the  potential  energy  operator  to  become  undiagonalized;  this  form  of  the 
Hamiltonian  is  called  the  non-adiabatic  representation.  The  coupling  angle  by  which  the 
adiabatic  representation  is  rotated  into  the  diabatic  representation  is  given  by  a  line 
integral  from  an  arbitrary  zero  to  the  configuration  in  question  through  the  NAD  terms. 
Non-adiabatic  chemistry  requires  a  quantum  chemistry  software  package  that  calculates 
NAD  terms.  Computational  results  from  two  packages,  Columbus  and  Brooklyn,  are 
compared  and  discussed. 

Separation  of  internal  dynamics  characterized  by  Jacobi  coordinates,  and  external 
dynamics  characterized  by  a  set  of  Euler  angles  and  the  center  of  mass  position,  requires  a 
transformation  from  Cartesian  coordinates,  employed  by  both  Columbus  and  Brooklyn,  to 
Jacobi  coordinates  required  for  subsequent  dynamical  calculations.  Previous  attempts  to 
solve  for  non-adiabatic  energy  surfaces  in  this  manner  have  failed  because  of  an 
ambiguity  in  selecting  the  correct  variable  for  describing  the  overall  rotation  of  the  B+H? 
system,  giving  answers  that  do  not  agree  with  specific  test  cases  for  which  the  coupling 
angle  is  known  via  simple  symmetry  arguments.  This  error,  which  lies  within  the  method 
of  converting  NAD  terms  from  one  coordinate  system  to  another,  is  discovered  and 
corrected.  By  way  of  this  correction,  correct  coupling  angles  are  determined,  and  non- 
adiabatic  energy  surfaces  are  calculated. 
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NON  -  ADI  AB  ATIC  ENERGY  SURFACES  OF  THE  B+H2  SYSTEM 


I.  Introduction 

Non-Adiabatic  Chemistry 

The  study  of  quantum  chemistry  often  takes  advantage  of  a  number  of 
approximations,  among  them  the  Bom-Oppenheimer  Approximation.  The  Born- 
Oppenheimer  Approximation  assumes  that  the  difference  in  mass  of  nuclei  and  electrons 
(the  smallest  ratio  is  about  1800:1  for  a  hydrogen  nucleus)  causes  them  to  have  vastly 
different  timescales,  such  that  the  dynamics,  and  thus  the  Hamiltonians,  of  the  two  groups 
can  be  separated.  Solution  of  the  electronic  Schrodinger  wave  equation  leads  to  a  set  of 
eigenvalues  and  eigenfunctions,  where  the  eigenvalues  serve  as  potential  energy  surfaces 
in  the  nuclear  Hamiltonian.  The  nuclear  Hamiltonian  then  takes  the  form  of  a  diagonal 
potential  energy  operator  and  a  kinetic  energy  operator  with  off-diagonal  terms.  This  form 
of  the  Hamiltonian  is  known  as  the  adiabatic  representation.  A  kinetic  energy  operator 
with  off-diagonal  elements  does  not  allow  for  solvable  differential  equations  for  the 
dynamics  of  the  system.  Fortunately  these  off  diagonal  elements  include  derivatives  of 
the  electronic  wave  functions  with  respect  to  the  nuclear  coordinates  (called  derivative 
coupling  terms ,  or  NAD  terms).  Because  the  Bom-Oppenheimer  Approximation  assumes 
different  electronic  and  nuclear  timescales  these  derivatives  can  be  approximated  to  zero, 
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and  both  operators  can  be  considered  diagonal  in  this  representation.  This  branch  of 
quantum  chemistry  is  called  adiabatic  chemistry. 

For  some  systems  the  dynamics  of  the  electrons  and  nuclei  are  not  so  cleanly 
separable.  We  may  still  use  the  Born-Oppenheimer  Approximation  to  separate  the 
nuclear  and  electronic  Hamiltonians;  nevertheless,  the  derivative  coupling  terms  are  not 
insignificant,  and  thus  the  Born-Oppenheimer  Approximation  cannot  be  used  to  assume  a 
diagonal  kinetic  energy  operator.  Since  a  diagonalized  kinetic  energy  operator  is 
necessary  for  creating  solvable  differential  equations,  additional  steps  must  be  taken  to 
manipulate  the  Hamiltonian  into  a  usable  form.  This  new  manipulated  representation  of 
the  Hamiltonian  that  diagonalizes  the  kinetic  energy  operator  (but  undiagonlizes  the 
potential  energy  operator)  is  called  the  non-adiabatic  representation,  and  consequently 
this  branch  of  quantum  chemistry  is  known  as  non-adiabatic  (or  diabatic)  chemistry. 

The  Scattering  Matrix 

In  non-adiabatic  chemistry,  diagonalization  of  the  kinetic  energy  operator  leads  to 
a  new,  undiagonalized  potential  energy  operator  and  new  associated  potential  energy 
surfaces.  The  goal  of  this  thesis  is  to  calculate  those  potential  energy  surfaces;  however, 
this  is  not  the  complete  picture.  These  surfaces  are  only  a  tool  with  which  to  construct 
the  scattering  matrix  which  describes  the  dynamics  of  the  system.  When  constructed,  the 
scattering  matrix  predicts  probabilities  of  reactions  occurring  at  the  molecular  level,  and 
thus  can  be  utilized  to  predict  how  chemical  reactions  will  proceed.  An  overview  of  the 
construction  of  the  scattering  matrix  appears  in  Appendix  A. 
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Relevance  to  Air  Force 


The  heavy  dependence  of  the  Air  Force  on  chemical  studies  for  materials,  fuels, 
meteorology,  etc.  makes  the  study  and  exploration  of  non-adiabatic  chemistry  a  necessary 
step  in  present  and  future  research.  For  example,  open-shell  collisions  such  as  N+CF  or 
O+No  that  affect  the  composition  of  the  upper  atmosphere  may  not  be  adequately 
described  by  approximating  the  derivative  coupling  terms  to  zero  in  the  adiabatic 
representation,  and  consequently  must  be  studied  in  the  non-adiabatic  representation. 

The  focus  of  this  study  is  another  system  of  interest.  Boron-doped  cryogenic 
hydrogen  (B+H2),  which  has  an  application  as  a  rocket  fuel  (Yarkony,  1999:i-2).  The 
lowest  energy  surface  (ground  state)  of  the  system  has  a  potential  well,  meaning  the 
system  can  assume  a  specific  configuration  and  will  not  be  able  to  change  without  outside 
energy.  This  well  corresponds  to  the  configuration  in  which  the  boron  atom  is  close  to 
the  center-line  of  the  hydrogen  molecule,  but  does  not  break  the  H-H  bond.  At  cryogenic 
temperatures  the  system  can  assume  this  configuration,  storing  energy  for  later  extraction 
when  the  fuel  is  burned.  This  well  and  the  associated  stored  energy  allow  B+FF  to  be 
classified  as  a  High  Energy  Density  Material  (HEDM).  B+H2  also  has  another  low- 
energy  surface,  close  to  but  higher  than  that  ground  state,  which  is  anti-bonding  (meaning 
there  is  no  potential  well).  Using  the  Bom-Oppenheimer  Approximation  to  neglect  the 
NAD  terms,  the  adiabatic  calculation  of  these  surfaces  assumes  that  the  bonding  surface 
and  the  anti-bonding  surface  are  not  associated;  that  is,  that  a  system  in  the  potential  well 
of  the  bonding  surface  will  not  couple  to  the  higher  anti-bonding  surface.  However,  if 
there  was  significant  coupling  between  the  surfaces,  the  B-UU  molecule  could  possibly 
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leave  the  well  and  dissociate,  shedding  its  stored  energy  prematurely,  rendering  it 
ineffective  as  a  high  energy  density  fuel.  The  degree  to  which  this  non-coupling 
assumption  is  true  is  not  easily  extracted  from  the  Hamiltonian  in  the  adiabatic  form. 
Rather,  that  information  lies  in  the  coupling  surface — the  collection  of  off-diagonal 
elements  of  the  non-adiabatic  potential  energy  operator.  Hence,  in  the  case  of  B+Hi,  one 
must  construct  the  non-adiabatic  representation  of  the  Hamiltonian  in  order  to  evaluate 
the  probability  of  coupling,  and  thus  suitability  as  a  high  energy  density  fuel. 

Recent  Work  and  Problems 

Dr.  David  R.  Yarkony,  Johns-Hopkins  University,  has  calculated  the  adiabatic 
surfaces  for  the  B-UU  system  as  well  as  the  derivative  coupling  terms  through  his  own 
software  package  called  Brooklyn.  Brooklyn  is  a  Graphical  Unitary  Group  Approach- 
(GUGA)  based  software  package  that  calculates  energy  levels  and  molecular  wave 
functions  on  the  Multireference  Configuration  Interaction  (MR-CI)  level  of  theory.  At 
the  time  of  publication  Brooklyn  was  not  available  for  public  use;  an  alternative  code, 
Columbus  developed  in  Columbus,  Ohio  originally  by  I.  Shavitt  et  al.  was  available  for 
public  use,  although  certain  portions  of  the  software  suite  are  still  under  development 
(Lischka,  undated-a).  Since  these  software  packages  handle  the  creation  of  molecular 
wave  functions  and  adiabatic  energy  surfaces,  they  can  be  seen  as  a  black  box;  it  is  not 
necessary  that  the  reader  understand  the  underlying  theory.  Nevertheless,  a  brief 
overview  of  methods  to  solve  the  many-body  problem  is  covered  in  Appendix  B. 

Since  Brooklyn  and  Columbus  are  designed  to  handle  a  wide  variety  of  molecules 
their  inputs  and  outputs  for  nuclear  positions  are  in  terms  of  Cartesian  coordinates.  For 
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most  systems  including  B-HT  we  would  prefer  to  develop  a  coordinate  system  that  allows 
for  the  separation  of  internal  and  external  dynamics  (and  thus  internal  and  external 
Hamiltonians).  In  the  case  of  the  three-body  system  the  internal  coordinates  are  the  Jacobi 
coordinates,  where  r  gives  the  bond  distance  of  Ho,  y  gives  the  tumbling  angle  of  boron 
with  respect  to  the  hydrogen  molecule,  and  R  gives  the  distance  between  boron  and  the 


Figure  1.  Jacobi  Coordinates 

center  of  mass  of  the  hydrogen  (see  Figure  1).  The  angle  y  can  be  measured  either  as  the 
angle  shown  in  Figure  1  or  its  supplement;  for  consistency  we  will  always  call  the  lesser 
of  the  two  angles  y.  Because  the  two  hydrogen  atoms  are  indistinguishable  it  is  not 
necessary  to  examine  the  entire  cycle  of  y  from  0  to  2n.  Throughout  the  calculations  we 
will  only  evaluate  y  from  0  to  tt/2  because  the  wave  function  will  be  symmetrical  in  all 
four  quadrants.  Another  advantage  of  using  Jacobi  coordinates  is  that  we  can  compare 
the  results  of  these  software  packages  with  the  work  of  Dr.  Millard  Alexander,  who 
calculated  some  non-adiabatic  potential  energy  surfaces  of  B-HT  in  Jacobi  coordinates  for 
r=  1.402,  the  equilibrium  bond  length  of  the  FT  molecule  (Alexander,  1993:6019). 
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Dr.  David  Weeks  took  upon  himself  the  job  of  transforming  Dr.  Yarkony’s  NAD 
terms  from  Brooklyn  into  their  Jacobi  counterparts.  These  results  did  not  match  those  of 
Dr.  Alexander.  Finding  the  correction  has  been  the  thrust  of  this  thesis,  since  we  cannot 
construct  the  proper  non-adiabatic  potential  energy  surfaces  without  knowing  the  correct 
NAD  terms.  To  find  and  correct  the  error,  this  study  was  split  into  two  explorations 
before  being  able  to  confidently  construct  the  non-adiabatic  potential  energy  surfaces: 
first,  determining  if  the  NAD  terms  produced  by  Brooklyn  were  correct  by  attempting  to 
reproduce  them  with  Columbus;  and  second,  determining  if  the  method  of  converting 
Cartesian  NAD  terms  to  their  Jacobi  counterparts  was  flawed.  The  final  result  is  that  the 
Brooklyn  NAD  terms  were  correct  and  it  was  the  conversion  process  that  contained  the 
error. 
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II.  Theory 


The  Hamiltonian 

As  stated  in  the  introduction,  the  goal  of  this  study  is  to  find  the  non-adiabatic 
potential  energy  surfaces,  specifically  for  B+FU.  Energy  surfaces  are  no  more  than  a 
collection  of  energy  eigenvalues  each  of  which  is  specific  to  a  particular  atomic 
configuration.  To  find  energy  eigenvalues  we  begin  at  the  time  independent  Schrodinger 
equation  (TISE)  (Liboff,  1998:72): 

(h-e)->P  =  0.  (1) 


The  form  expressed  in  equation  (1)  is  the  most  basic  and  abstract  way  to  view  the  action 
of  the  Hamiltonian  operator.  To  be  of  any  use,  we  must  consider  it  in  some  coordinate 
system  and  apply  it  to  a  molecule,  at  least  generally.  The  Hamiltonian  operator  for 
polyatomic  molecules  in  the  coordinate  representation  is  written  as: 


M, 


“ZV-Z- 

i  fY 


R, 


P\Ra-Rt 


(2) 


where  i  and  j  index  electrons,  a  and  P  index  the  nuclei,  r  is  the  electronic  coordinates  and 
R  is  the  nuclear  coordinates  (Szabo  1996:41).  Unless  otherwise  noted  we  are  using 
atomic  units  (au).  The  terms  are: 
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r  -  R, 


Z 
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z  z 


Nuclear  Kinetic  Energy 


Electron  Kinetic  Energy 


Nucleus-Electron  Interaction 


Electron-Electron  Interaction 


Nucleus-Nucleus  Interaction. 


The  Born-Oppenheimer  Approximation 

The  Born-Oppenheimer  approximation  assumes  that  the  nuclei  move  far  slower 
than  the  electrons  (Bransden,  1984:386-88).  Thus  to  the  electrons,  the  nuclei  can  be  seen 
as  standing  still,  creating  a  constant  (time-independent)  potential  field.  To  the  nuclei  the 
electrons  move  so  fast  that  their  average  can  be  seen  as  a  constant  potential  field.  For 
electron  dynamics  this  approximation  allows  us  to  assume  that  the  first  term  of  the 
Hamiltonian  (nuclear  kinetic  energy)  is  equal  to  zero,  and  the  last  term  (potential  due  to 
intemuclear  Coulombic  forces)  is  equal  to  a  constant.  The  effect  of  adding  a  constant  to 
the  Hamiltonian  is  only  to  add  a  constant  to  the  eigenvalues;  since  we  can  arbitrarily 
choose  the  zero  of  a  potential  field  we  can  neglect  this  constant  term.  The  remaining 
terms  form  the  electronic  Hamiltonian: 


i,a  ri  ^ a  i,j 


(3) 


r.  —  r  ■ 

'  J I 


which  has  a  complete  set  of  orthogonal  electronic  eigenfunctions,  cp;(R;  r).  Note  that, 
although  the  nuclear  terms  are  left  out  of  the  Hamiltonian,  the  nuclear  configuration  is  not 
to  be  ignored.  The  functions  cpi(R;  r)  depend  parametrically  upon  the  nuclear  coordinates 
(R)  as  much  as  they  depend  functionally  on  the  electronic  coordinates  (r).  For  simplicity 
I  will  use  abstract  notation  until  necessary  to  operate  in  the  coordinate  representation. 

The  eigenfunctions,  (pi,  of  the  electronic  Hamiltonian  form  a  complete  and 
orthogonal  set,  so  that  any  function,  specifically  the  wave  function,  \\i,  from  equation  (1) 
can  be  written  as  a  superposition  electronic  wave  functions: 

W-Zk-k)  <4> 

i 

where  the  value  of  the  qth  component,  Fq  ,  can  be  found  by: 

(n  |  'p) = I  r,  (n  \<p, ) = fa,, =Ft,  <5) 

i 

where  the  orthogonally  of  the  electronic  states  (pi  has  been  used.  Since  the  total 
Hamiltonian  is  a  product  of  the  nuclear  and  electronic  Hamiltonians  (under  the  Born- 
Oppenheimer  Approximation),  the  eigenfunctions  of  the  full  Hamiltonian  can  be  written 
as  a  product  of  eigenfunctions  of  those  Hamiltonians.  Thus  if  |  T*)  from  equation  (4)  is  an 

eigenfunction  of  the  total  Hamiltonian  and  |  (pt )  are  the  eigenfunctions  of  the  electronic 

Hamiltonian,  it  follows  that  the  F,’s  are  the  eigenfunctions  of  the  nuclear  Hamiltonian, 
and  as  such  are  functions  of  the  nuclear  coordinates  (this  will  be  key  to  the  approximation 
shortly). 
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Since  for  the  electronic  Hamiltonian 

Z  (n  |h,  k ) = Z  (n  K  k) = = £, .  (6) 

i  i 

it  follows  that  for  the  total  Hamiltonian, 

HT=He+HN,  (7) 

(where  H N  is  the  nuclear  Hamiltonian)  we  have  the  following  derivation.  We  start  by 
integrating  the  TISE  in  equation  (1)  with  respect  to  electronic  coordinates: 

(nN'f'HnN'p)-  <8) 

We  then  substitute  for  the  total  Hamiltonian  with  equation  (7)  and  substitute  for  v| /  with 
equation  (4), 

Z  k  |h. + «»|  fm) = Z(n  N  F&<)  (»> 

i  i 

which  can  then  be  split: 

Zk  |kkk+Zk  lAk^-MZk  |r«>-  ao 

i  i  i 

And,  taking  advantage  of  (6)  we  have 

E,F,+'L{‘p,'?iAFt‘p,)=EF,  do 

i 

or 

ZkN  F,<P,)  +  (F,-E)-F,  =0.  (12) 

i 

We  do  not  yet  know  how  to  evaluate  the  term  ^  ((pq  \H  N  |  Fi(pi ) .  In  order  to  evaluate  it  we 

i 

must  consider  it  in  the  coordinate  representation.  Recall  that 
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(13) 


2  M 

^  a  1¥1  a 


thus,  when  operating  on  the  full  molecular  wave  function  v| /, 


i 


(14) 

where  we  have  merely  used  the  product  rule  for  derivatives.  Notice  the  pair  of  terms 
(vAy  F  (/?))•  (pi  (R;  r)  and  (Vs$>.  (R:  r)).  These  terms  take  derivatives  of  the  electronic  wave 

functions  with  respect  to  nuclear  coordinates.  These  terms,  when  integrated  with  respect 
to  electronic  coordinates  become  the  derivative  coupling  terms  set  forth  in  the 
introduction.  The  Born-Oppenheimer  Approximation  allows  us  to  assume  that  derivatives 
of  the  electronic  wave  function  with  respect  to  the  nuclear  coordinates  are  negligible;  that 
is,  Va(pq  and  Va2cpq  are  zero,  so 

(r\H„\V)  =  H„ZFl(R)MKr)  =  2Z^{{v/FtR))-pl{R-,r)}.  (15) 

^  a,i  At  a 

This  leads  to  another  uncoupled  TISE.  If  we  now  substitute  the  result  of  (15)  into  (12) 
we  get: 

|Z-^-V»A+(£,-Af,=0  (16) 

L  a  1V1  a 


or,  rearranging  the  terms  into  an  eigenvalue  equation, 


-x— 

y2  a  M  a 


-vs-+£9 


■  f(r)=ef(r). 


(17) 
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Equation  (17)  is  the  nuclear  TISE,  with  Fq(R)  as  its  eigenfunctions.  Eq(R),  the 
eigenvalues  for  the  electronic  wave  equation,  now  become  the  effective  potential  for  the 
nuclear  wave  equation.  The  Hamiltonian  is  broken  into  two  operators,  kinetic  and 
potential  energy,  which  can  be  represented  as  matrices  operating  on  the  nuclear  wave 
function.  Since  the  TISE  is  an  eigenvalue  equation,  we  can  refer  to  the  nuclear  wave 
functions  as  its  eigenvectors  or  eigenfunctions. 

Using  the  Bom-Oppenheimer  approximation,  forming  matrix  operators  out  of  the 
Hamiltonian  is  simple.  Before  we  do,  however,  we  must  discuss  truncation  of  the  energy 
surfaces. 

Basis  and  Symmetry 

For  the  B+H2  system,  there  are  an  infinite  number  of  molecular  orbitals,  each  of 
which  has  an  associated  energy  surface,  that  could  be  taken  into  account  (Szabo,  1996:55- 
57).  These  molecular  orbitals  serve  as  basis  vectors  for  forming  the  nuclear  wave 
functions.  It  is  of  course  impossible  to  build  the  infinite  matrix  associated  with  such  a 
basis.  Even  a  small  number  of  these  surfaces  will  begin  to  present  numerical  difficulties 
(Szabo,  1996:58).  In  order  to  keep  calculations  to  a  minimum  but  still  arrive  at 
worthwhile  results,  we  will  only  consider  the  molecule  in  its  electronic  ground  state.  At 
the  very  low  temperatures  at  which  this  system  is  of  interest  as  a  cryogenic  fuel,  these  are 
the  only  energy  levels  of  interest. 

A  peculiarity  of  this  molecule  is  that  the  valence  electron  of  boron  has  three 
choices  of  atomic  orbital  for  its  ground  state  owing  to  the  three-fold  degeneracy  of  the  p- 
orbital  (Alexander,  1996:6015).  At  the  asymptotic  limit,  when  boron  is  very  far  from 
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hydrogen,  the  eigenvalues  of  these  three  choices  of  orbital  are  degenerate;  but  as  FF  and 
B  come  together  the  energy  surfaces  associated  with  these  three  ground  states  begin  to 
separate,  but  are  still  closer  together  than  higher  states.  For  this  reason  we  will  only 
consider  the  three  energy  surfaces  that  arise  from  the  system’s  occupying  each  of  three 
molecular  orbitals,  each  of  which  is  composed  primarily  of  only  one  of  boron’s  atomic  p- 
orbitals  for  our  energy  levels  of  interest.  Thus  all  the  nuclear  wave  functions  we  consider 
will  be  linear  combinations  of  these  three  molecular  orbitals. 

Special  note  should  be  taken  that  since  B+FF  has  only  three  atoms,  it  can  always 
be  placed  in  a  plane.  Thus,  except  for  the  linear  and  perpendicular  configurations,  the 
system  has  the  highest  symmetry  of  Cs,  meaning  the  only  non-identity  symmetry  it  has  is 
reflection  through  the  plane  of  the  molecule  (Bishop,  1993:37).  The  molecular  orbitals 
can  thus  be  symmetric  across  the  plane,  which  are  given  a  symmetry  label  A’  or 
antisymmetric,  which  are  given  a  symmetry  label  A”.  Figure  2  shows  an  example  of 
symmetric  and  antisymmetric  orbitals  with  respect  to  the  x-y  plane.  As  it  turns  out  the 
three  ground  states  for  the  system  include  one  wave  function  of  A’  ’  symmetry  labeled 


Figure  2.  Symmetric  vs.  antisymmetric  p-orbitals  when  reflected  through  the  x-y  plane 
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1A”,  and  two  of  A’  symmetry  labeled  1A’  and  2 A’  (Alexander,  1993:6018).  Each 
molecular  orbital  has  an  associated  potential  energy  surface  which  carries  the  same  label. 

If  we  consider  three  possible  molecular  orbitals  then  equation  (17)  in  matrix  form 
becomes: 


Px  1  v*2 

2  ^  M 

0 

0 

“1A' 

0 

0  1 

0 

‘y  1  v„3 

2 

0 

+ 

0 

2  A' 

0 

a 

1  1 

0 

0 

1A"J 

0 

0 

y  v„2 

2  rMa  \ 

(18) 

(Weeks,  2004). 


Note  that  the  number  of  molecular  orbitals  we  wish  to  consider  determines  the 
dimensions  of  the  operator  matrix.  Hence  we  can  think  of  doing  algebra  in  a  three-space, 
where  the  three  molecular  orbitals  serve  as  basis  vectors.  Since  the  potential  matrix  has 
no  off-diagonal  elements,  this  is  the  adiabatic  representation.  This  condition  assures  that 
a  nuclear  wave  function  on  one  potential  energy  surface  does  not  interfere  with  the 
function  on  another.  Since  the  Bom-Oppenheimer  Approximation  allows  that  the  off- 
diagonal  elements  of  the  kinetic  energy  operator  be  assumed  negligible,  it  is  diagonal  as 
well  and  we  can  extract  solvable  equations  from  this  eigenvalue  equation  where  nuclear 
dynamics  on  each  of  the  1A’,  2A’,  and  1A”  surfaces  is  uncoupled. 


Refining  the  Approximation 

We  will  now  consider  the  situations  in  which  the  Born-Oppenheimer 
Approximation  cannot  be  used  to  go  from  equation  (14)  to  equation  (15).  That  is,  the 
derivatives  of  the  electronic  wave  functions  with  respect  to  the  nuclear  positions  or  NAD 


14 


terms  are  significantly  greater  than  zero,  and  must  be  included.  Recall  that  these  terms 
are  the  off-diagonal  terms  of  the  kinetic  energy  matrix  in  equation  (18).  When  these 
terms  are  included,  the  potential  energy  matrix  remains  diagonalized,  but  the  kinetic 
energy  matrix  is  not.  The  result,  from  equation  (14),  is 


-Z— V/ 

^  a  lvla 

0 

0 


-S— V/ 

^  a  iV1a 

o 


0 

0 


-S— V«2 

^  a  iV1a 


+ 


ZrrWkfl  *va.+v,^]  Z^Pr^lv/^'vA'+vA'2^]  Zkpr$k$  *vA,+vA,vJ 


a  lvla 


a  1Y1a 


a  lvla 


Zxrf^’kfl  *va+va2^]  Z^Pr^lvA^'vA+vA'Vj  *VA,+VA,Vj 


a  iV1a 


a  lvla 


a  lvla 


ZxrK^k^?  •  +v,2*]  z^-f rf^k&  •  v« +V^J  Zjkk4v^#v* + vA2^] 


a  1Y1a 


■'M 

a  ,VIa 


‘M 

a  lyla 


+ 


X 

0 

o'1 

pT 

X 

0 

e2 

0 

A 

=E 

F 

0 

0 

A 

Ia 

A 

•(19) 

Upon  inspection  it  can  be  seen  that  this  is  the  same  as  equation  (18),  but  to  the  kinetic 
energy  operator  has  been  added  a  3x3  matrix  of  derivative  coupling  operators.  To  further 
simplify  this  representation,  we  recognize  that 

P„{R)=jdrf>;(R-,r)-V,9l(R-,r)  (20) 

is  just  a  constant  dotted  into  the  del  operator  and 

Qf(R)=jdrf>;{R-,r)^2,f>l{R-,r)  (21) 
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is  also  a  constant.  Note  that  the  differential  operator  is  with  respect  to  nuclear 
coordinates  but  the  integral  is  with  respect  to  electronic  coordinates.  These  are  the 
derivative  coupling  terms.  Hence  we  can  write  the  Hamiltonian  as 


J_y  1  . 

2  rMa 


0 

0 


-X— V/ 

^  a  1V1  a 


o 

0 


-E— 

O  A/I 


2  rMa 


+ 


EvJ-k-V.+Qj  EtM^^  +  qJ 

M  M 

a  1V1  a  a  1V1  a 

vfi+e21]  E^k  vff  +  G22] 


a  M  a 


a  M  a 


E{  0  0 

0  e2  0 
0  0  e2 


0 

0 


E  A/r  [^*33  R  +  033  ] 


a  M  a 


+  ■ 


F, 

E 2 
F, 


=  E 


F 2 

F, 


(22) 


Ordinarily  the  13,  23,  31,  and  32  elements  of  the  second  matrix  would  be  non-zero  as 
well;  however,  the  1A’  and  2 A’  molecular  orbitals,  which  are  primarily  composed  of  the 
atomic  p-orbitals  of  boron  that  are  coplanar  with  the  system  tend  to  mix  with  each  other, 
but  not  with  the  1A”  molecular  orbital,  which  is  composed  primarily  of  the  atomic  p- 
orbital  orthogonal  to  the  plane.  Thus  we  save  CPU  time  by  recognizing  that  those 
elements,  which  represent  mixing  of  the  third  with  either  of  the  first  two,  are  nearly  zero 
(Alexander,  1993:6019). 


The  Diabatic  Transformation 

If  we  were  satisfied  with  the  adiabatic  representation  we  could  stop  here; 
unfortunately  without  employing  the  Bom-Oppenheimer  Approximation  this 


16 


representation  does  not  provide  a  set  of  differential  equations  which  are  solvable.  If, 
however,  the  kinetic  energy  matrices  are  considered  in  a  different  basis  in  which  they  are 
diagonalized  we  will  form  three  new  diabatic  potential  energy  surfaces  and  a  set  of 
diabatic  equations  that  are  tractable. 

In  order  to  understand  how  to  diagonalize  the  kinetic  energy  operator,  consider 
that  the  three  molecular  orbitals  1A’,  1A”,  and  2A’  are  the  basis  vectors  of  a  wave 
function  three-space.  Consider  a  nuclear  wave  function,  F,  to  be  a  vector  in  this  three- 
space,  as  depicted  in  Figure  3.  One  could  imagine  splitting  this  vector  into  three 
components,  one  along  each  of  the  axes.  This  depiction  corresponds  to  constructing  the 
wave  function  from  a  linear  combination  of  these  molecular  orbitals,  or  finding  parts  of 


Figure  3.  A  Nuclear  Wave  Vector  Represented  in  Wave  Function  Three-Space 

the  total  wave  function  on  all  three  potential  energy  surfaces.  A  linear  operator,  O,  that 
acts  on  F  will  act  on  the  components  of  F  individually,  as  one  might  expect  simple  in 
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linear  algebra.  If  6  is  diagonal,  then  the  result  of  O(Fi)  is  also  along  component  i  (where  i 
is  1A’,  1A”,  or  2A’).  If  O  is  not  diagonal,  but  there  exist  three  eigenvectors  of  O  which 
form  an  alternative  basis  for  the  three-space,  then  O  is  diagonalizable  in  that  basis 
(Hoffman,  1971:185).  The  diagonalization  takes  the  form  of  pre-  and  post-multiplication 
by  a  rotation  operator,  T  (created  with  knowledge  of  the  inner  products  of  the  original 
basis  vectors  and  the  eigenvectors  of  O)  such  that 

dD=TOTx  (23) 

where  0D  is  the  representation  of  6  in  the  new  basis,  and  is  diagonal.  The  basis  vectors 
of  this  new  basis  are  the  eigenvectors  of  O  used  to  create  the  rotation  operator.  They  will 
be  some  mixture  of  the  original  basis  vectors,  as  illustrated  in  Figure  4.  This  mixing 


Figure  4.  Rotation  of  Potential  Energy  Surface  Basis 

depicts  the  change  from  molecular  orbitals  and  potential  energy  surfaces  in  the  adiabatic 
representation  to  molecular  orbitals  and  potential  energy  surfaces  in  the  non-adiabatic 
representation.  The  1A’,  2A’  and  1A”  vectors  are  eigenfunctions  of  the  adiabatic 
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potential  energy  operator,  and  the  1AD\  2AD’  and  1AD”  vectors  are  the  eigenfunctions  of 
the  diabatic  kinetic  energy  operator. 

From  equation  (22)  we  know  that  our  kinetic  energy  operator  which  we  shall  call 
K  is  not  diagonal.  Creating  the  rotation  operator  matrix  to  diagonalize  K  would  normally 
require  three  angles  of  rotation  that  mix  the  old  molecular  orbitals  into  new  ones.  But 
since  the  1A”  orbital  does  not  mix  with  the  others,  we  need  only  worry  about  mixing  the 
two  A’  orbitals  which  requires  only  a  single  angle  0(R)  (Koppel  2002:7).  This  rotation, 
applied  the  matrix-form  Hamiltonian,  will  reduce  the  P(R)  and  Q(R)  terms  to  zero.  The 
rotation  matrix  takes  the  form 


such  that 


T  = 


cos  {®{r)) 
-  sin(0(i?)) 
0 


sin(0(7?)) 

cos(0(/?)) 

0 


0 

0 

1 


K  =  TKT  X\U  =  TUT1 


(24) 


(25) 


where  U  is  the  potential  energy  operator  and  the  primed  matrices  are  the  kinetic  and 
potential  energy  operators  in  the  diabatic  representation  In  this  representation, 
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(26) 
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(27) 


E\ , ,  E12 ,  and  are  the  diabatic  (or  non-adiabatic)  potential  energy  surfaces  upon 
which  we  can  propagate  the  wave  functions;  and  Eu ,  though  not  a  proper  potential 
surface,  is  the  coupling  value  between  the  first  two  surfaces.  [K’+U’]F=EF  is  now  a 
solvable  system  of  equations,  and  furthermore  we  have  a  measure  of  coupling  between 
the  potential  energy  surfaces. 

The  task  now  turns  to  identifying  0(R),  which  will  depend  on  P(R)  and  Q(R).  As 
an  approximation,  let  us  assume  that  Q(R),  which  is  a  double  derivative  of  electronic 
eigenfunctions  with  respect  to  nuclear  coordinates,  is  negligible  compared  to  P(R)  which 
is  only  a  single  derivative  of  the  same.  Recall  from  equation  (20)  that  Py  has  the  form 

PII(R)  =  \dr<pl(R-r)VR<pl(R-,r)  (20) 

or,  abstractly, 

Pij={(Pi\VR(Pj).  (28) 

In  what  Horst  Koppel  calls  the  “off-diagonal  analogue  of  the  Hellmann-Feynman 
theorem”  we  can  manipulate  the  NAD  term  into  a  more  useable  formula  (2002:5).  Let  us 
start  again  with  the  TISE: 

«k)=p  h)-  <29) 

We  then  take  the  gradient  with  respect  to  the  nuclear  coordinates: 

v(h|^})=v(£.|^))  (30) 
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which  by  the  product  rule  becomes 


V^|^)  +  ^|V^)  =  V^|^)  +  ^|V^).  (31) 

Now  integrate  this  equation  over  electronic  coordinates: 

(^.  |VH  |  ^.  )  +  (?>.  |  V  <p.  )  =  {(Pi  |  VEj  |  <Pj)  +  (<pt  \Ej\V  (Pj ) .  (32) 


Both  energy  terms  on  the  right-hand  side  of  the  equation  are  constants  and  can  be 
removed  from  the  integration,  while  the  Hamiltonian  in  the  second  term  on  the  left-hand 
side  of  equation  (32)  can  act  on  the  bra  rather  than  the  ket: 

{<Pi  Fh\<Pj)  +  Ei {<Pi\v <Pj )  =  VEj (<P>  \<Pj)  +  ej M V<P}).  (33) 

We  now  have  an  equation  which  has  the  NAD  term  in  it.  Elimination  of  the  first  term  on 
the  right-hand  side  by  orthogonality  and  simple  algebraic  manipulation  reveal  a  new 
formula  for  the  NAD  term: 


pa=\<Pi  V<Pj)  = 


(Pi  VH  <p 


Ej-Ei 


(34) 


For  an  overview  of  how  Columbus  calculates  Py,  see  Appendix  C.  This  form  of  the 
derivative  coupling  term  brings  to  light  one  situation  where  the  Bom-Oppenheimer 
approximation  cannot  be  applied  to  the  adiabatic  representation.  In  the  denominator  of 
the  term  is  the  difference  between  the  two  energy  surfaces  i  and  j.  As  they  become  close, 
the  derivative  coupling  term  becomes  large  and  cannot  be  neglected.  Moreover,  at  the 
point  where  the  surfaces  cross  the  term  becomes  a  singularity,  creating  havoc  with 
adiabatic  predictions. 

Advantage  may  be  taken  upon  examining  the  symmetry  of  the  electronic  wave 
functions  in  the  NAD  term.  Since  V  is  a  symmetric  operator,  if  one  of  the  wave  functions 
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integrated  in  Py  is  symmetric  but  the  other  is  antisymmetric,  the  entire  integral  is 
antisymmetric  over  all  space,  and  is  zero.  This  is  why  no  coupling  is  possible  (at  this 
level  of  theory)  between  the  A’  and  A”  states,  and  thus  why  when  we  consider  the 
derivative  coupling  terms,  we  shall  only  consider  the  terms  that  couple  1A’  and  2A\ 

The  new  P’,j  of  the  diabatic  representation  will  be  related  to  the  adiabatic  P,,  by  a 
gauge  transformation  (Koppel,  2002:7): 

P;i(R)=Pll(R)+V(0(R)).  (35) 

We  want  to  set  this  to  zero,  so 

V(q(r))  =  -Pij(r)=Pji(r)  .  (36) 

Solving  for  0(R)  simply  requires  integration  of  this  equation: 

R 

®{R)  =  -\dR'Pij{R')+Q{R0),  (37) 

o 

where  0(Ro)  represents  an  arbitrary  point  at  which  to  begin  the  line  integral. 
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III.  Calculation,  Results  and  Discussion 


Columbus  and  Brooklyn 

The  software  suite  of  choice  for  this  venture  was  Columbus.  Though  operation 
required  a  delicate  touch,  input  scripts  were  cryptic,  and  documentation  was  scarce,  in  the 
end  Columbus  was  the  only  package  available  to  the  public  that  could  compute  derivative 
coupling  terms  necessary  for  a  diabatic  transformation.  When  computing  energy 
eigenvalues  and  derivative  coupling  terms,  Columbus  implements  multireference 
configuration-interaction  (MR-CI)  (Lischka,  2004).  In  addition  to  Columbus ,  it  should  be 
noted  that  we  also  used  results  from  Brooklyn ,  proprietary  code  of  Dr.  Yarkony.  Since  Dr. 
Yarkony  has  aided  in  writing  Columbus  it  is  safe  to  assume  that  the  functions  that  the  two 
packages  have  in  common  are  very  similar. 

Columbus  itself  is  not  a  single  program,  but  a  suite  of  different  programs 
communicating  via  files.  This  modular  nature  allows  the  user  to  utilize  some  parts  of 
Columbus  without  having  to  bother  with  others.  In  order  to  get  derivative  coupling  terms 
as  we  needed,  Columbus  must  summon  up  a  dozen  or  more  individual  programs  to  create 
bases,  construct  molecular  wave  functions,  make  distinct  row  tables,  etc.  Each  one  of 
these  programs  individually  requires  a  number  of  input  files  which  are  not  easily 
constructed  by  hand.  For  this  reason  the  maintainers  of  Columbus  have  created  a  script 
file  COLINP.X  which  guides  the  user  through  choices  of  input  and  intended  output  and 
creates  the  necessary  files.  Like  the  rest  of  Columbus,  COLINP.X  is  a  work  in  progress 
and  still  requires  some  manual  tweaking  of  the  files  after  they  are  created.  It  also  creates 
a  file  input  for  RUNC.X,  another  script  file  that  will  call  the  necessary  programs  within 
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Columbus  that  are  needed  so  the  user  does  not  have  to  load  each  program  individually. 

More  information  about  the  operation  and  capabilities  of  Columbus  can  be  found 
at  the  Columbus  web  site  (Lischka,  undated-b). 

Coordinates 

Because  Columbus  and  Brooklyn  work  in  Cartesian  coordinates  it  was  necessary 
to  convert  our  Jacobi  coordinates  into  Cartesians  for  entry  into  the  software  packages. 
Since  each  of  the  three  atoms  must  have  its  position  fixed  in  space,  we  must  translate 
three  Jacobi  coordinates  into  nine  Cartesian  coordinates.  Obviously  the  molecule  can  be 
restricted  to  the  z=0  plane,  eliminating  the  need  for  three  extraneous  coordinates.  Within 
the  x-y  plane  the  three  Jacobi  coordinates  specify  how  the  atoms  are  oriented  with  respect 
to  one  another,  but  there  is  still  no  criterion  for  how  the  molecule  should  be  oriented 


y 


Figure  5.  Translating  Jacobi  Coordinates  into  Cartesian  Coordinates 
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within  that  plane.  At  the  request  of  Dr.  Yarkony,  a  decision  was  made  so  that  the  least 


number  of  Cartesian  coordinates,  three,  were  needed  to  specify  the  geometry  (Weeks, 
2004).  This  necessitated  the  pinning  of  boron  to  the  origin,  pinning  one  hydrogen  to  the 
x-axis  but  allowing  it  to  slide  along  it,  and  allowing  the  second  hydrogen  atom  to  be 
oriented  freely  in  the  plane  (see  Figure  5).  This  restriction  allowed  the  geometry  to  be 
specified  with  three  Cartesian  coordinates:  xHi,  xH2,  and  yH2-  Given  Jacobi  coordinates  r, 
R,  and  y,  the  Cartesian  coordinates  can  be  calculated  by 


where  the  sign  in  equation  (40)  is  equal  to  the  sign  of  the  dot  product  of  the  vector  along  r 
and  the  x  axis. 

Calculating  the  Adiabatic  Potential  Energy  Surfaces 

There  are  three  surfaces  of  interest  corresponding  to  the  three  molecular  orbitals 
1A’,  2A’  and  1A”  as  mentioned  previously.  A  slice  of  these  three  surfaces  as  calculated 
by  Columbus  is  shown  in  Figure  6.  Notice  that  as  R  goes  to  infinity,  that  is  as  the  boron 
atom  and  the  hydrogen  molecule  become  separated,  the  three  energy  surfaces  tend  to 
become  degenerate.  The  left  side  of  the  graph  shows  the  large  potential  energy  barrier  to 
becoming  BH2.  Figure  7  shows  an  expanded  view  of  the  flat  region  of  the  surfaces. 
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Notice  that  the  1A’  state  does  indeed  have  a  small  well  for  bonding,  whereas  the  2A’  state 
is  anti-bonding.  1A”  is  also  bonding,  but  cannot  be  coupled  to  the  others  because  of 
symmetry.  Both  of  these  graphs  show  the  slice  of  the  energy  surface  taken  at  y  =  n/4. 

The  deepest  part  of  the  well  actually  occurs  at  y  =  ji/2. 

Let  us  also  look  at  the  surface  as  sliced  through  y.  Figure  8,  Figure  9,  and  Figure 
10  all  show  slices  of  the  energy  surfaces  at  R=5.0  and  r=  1 .402  bohr  (equilibrium  distance 
for  hydrogen  molecule)  but  calculated  in  different  manners.  Figure  8  and  Figure  9  both 
show  the  same  slice  of  the  energy  surfaces  as  calculated  by  Columbus,  but  there  is  a 
difference  in  input.  Columbus  builds  molecular  orbitals  out  of  a  linear  combination  of 
atomic  orbitals.  In  theory  there  are  an  infinite  number  of  atomic  orbitals  available,  but  in 
practice  we  must  choose  a  finite  number  to  serve  as  a  basis  set  for  constructing  the 


Figure  8.  Energy  Surface  Slice  for  R=5.0  r=  1.402 
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Figure  9.  Energy  Surface  Slice  for  R=5.0  r=1.402,  Truncated  Basis 

molecular  orbitals  (this  is  not  to  be  confused  with  the  basis  of  molecular  orbitals  which 
form  a  basis  of  the  wave  function  three-space).  Figure  8  shows  the  potential  energy 
surfaces  as  calculated  with  the  triple  zeta  (PVTZ)  basis  set;  in  the  case  of  B+FE  this  set 
consists  of  58  basis  orbitals,  40  of  A’  symmetry  and  18  of  A”  symmetry.  Figure  9  shows 
the  same  potential  energy  surface  slice  as  calculated  with  the  double  zeta  (PVDZ)  basis 
set;  in  the  case  of  B+FE  this  set  consists  of  24  basis  orbitals,  18  of  A’  symmetry  and  6  of 
A”  symmetry.  Notable  differences  include  a  difference  in  energy  level  and  a  kink  in  the 
PVDZ  energy  curves.  The  PVTZ  energy  surfaces  start  below  the  PVDZ  set  by  at  about 
0.02  atomic  units.  The  variation  principle  states  that  the  energy  calculated  using  the  true 
Hamiltonian  is  always  an  upper  bound  to  the  true  energy,  suggesting  that  the  triple  zeta 
calculations  are  more  accurate.  Figure  10  shows  data  collected  from  Dr.  Yarkony’s 
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calculations  using 
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Figure  10.  Energy  Surface  Slice  for  R=5.0  r=1.402,  Brooklyn  Code 

Brooklyn.  As  best  we  can  tell  by  examining  his  input  files  he  is  using  83  atomic  orbitals 
in  his  basis  set.  This  could  suggest  mixing  PVTZ  and  PVDZ  bases  or  some  other  basis 
set.  Despite  the  larger  basis  set,  Dr.  Yarkony’s  energy  surfaces  seem  to  come  short  of 
Columbus’  PVTZ  calculations  by  about  0.005  atomic  units.  Although  either  of  these 
calculations  could  be  considered  sufficient,  it  may  suggest  that  the  method  Columbus  uses 
to  calculate  adiabatic  energy  surfaces  is  more  efficient  or  more  in-depth  than  that  of 
Brooklyn.  Figure  1 1  shows  surface  and  contour  plots  of  the  three  adiabatic  potential 
energy  surfaces  of  interest  as  calculated  with  Brooklyn  for  the  equilibrium  position  of 
hydrogen,  r=1.402,  for  the  range  R=5  to  10  and  y=0  to  jt. 
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Figure  11.  Surface  and  Contour  Plots  of  the  1A',  2A',  and  1A"  Surfaces 

Calculating  the  Derivative  Coupling  Terms 

Another  function  of  Columbus  or  Brooklyn  which  is  critical  to  calculating  the 
non-adiabatic  potential  energy  surfaces  is  calculation  of  the  derivative  coupling  terms. 
When  either  program  calculates  wave  functions,  it  only  does  so  to  within  a  phase  factor 
(Weeks,  2004).  When  the  electronic  wave  functions  are  integrated  to  get  the  NAD  terms, 
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this  phase  factor  is  integrated  out  for  diagonal  terms: 


P,  =  =  (41) 


But  for  off-diagonal  terms  where  0;  does  not  equal  0j  the  phase  factor  does  not  disappear. 
Fortunately  for  real  wave  functions  the  phase  factor  is  always  real,  so  the  phase  factor  just 
becomes  an  arbitrary  1  or  -1  in  front  of  every  calculated  NAD  term.  Unfortunately  there 
is  no  easy  fix  to  this  phase  factor,  but  neither  can  it  be  ignored.  The  only  way  to  correct 
for  it  is  to  view  the  NAD  terms  in  sequence  across  the  energy  surface  and  invert  terms 
that  seem  out  of  place.  The  two  guiding  facts  for  selecting  those  terms  in  need  of 
inversion  are  that  the  NAD  terms  should  vary  slowly  across  all  parts  of  the  surface  and 
that  the  phase  factor  is  the  same  for  all  NAD  terms  at  a  given  nuclear  configuration. 
Figure  12  illustrates  a  series  of  NAD  terms  calculated  by  Columbus.  The  label  d/dn  is 


shorthand  for  the  derivative  coupling  term  with  respect  to  the  coordinate  n:  (  cpj 


The  slice  is  the  same  as  previously  examined:  R=5.0  and  r=1.402  bohr.  Note  that  for  y= 
47i/20,  5tt/20,  and  6tt/20  the  all  of  the  terms  take  a  sharp  dive  across  the  axis.  If  it  is 
assumed  that  these  three  points  have  been  calculated  with  a  -1  phase  shift,  we  can  correct 
them  with  a  simple  sign  adjustment.  Figure  13  shows  the  same  NAD  terms  with  the 
phase  correction  made.  Columbus  calculated  the  points  marked,  whereas  the  lines  are 
interpolations  on  those  points.  We  could  have  chosen  the  exact  opposite  situation  and 
assumed  that  those  three  points  were  correct  and  all  the  others  were  flipped.  The  result 
would  be  the  same,  only  inverted  about  the  axis.  The  total  phase  factor  is  irrelevant  so 
long  as  it  is  the  same  for  all  measurements.  We  cannot,  however,  choose  to  correct  some 
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Figure  12.  Cartesian  NAD  Terms  from  Columbus  with  Uncorrected  Phase 
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Figure  13.  Cartesian  NAD  terms  from  Columbus  with  Corrected  Phase 
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of  the  NAD  terms  at  a  given  nuclear  configuration  and  not  others.  For  instance  in  Figure 
13  it  would  be  tempting  to  correct  the  d/dxhi  and  d/dyhi  terms  alone  at  y  =  71/IO.  This  is 
forbidden  without  flipping  the  other  four  terms  as  well,  and  by  examining  the  d/dxb  and 
d/dxh2  terms  we  can  judge  that  inversion  would  not  be  prudent.  Figure  14  shows  the 
NAD  terms  for  the  same  slice  of  potential  energy  surface  as  calculated  by  Dr.  Yarkony  in 
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Figure  14.  Cartesian  NAD  Terms  from  Brooklyn  with  Corrected  Phase 


Brooklyn  after  a  corrected  phase  factor  has  been  applied.  Unlike  for  the  potential  energy 
surface  calculations,  the  NAD  term  outputs  from  Columbus  and  Brooklyn  are  not  close  at 
all.  Not  only  do  they  have  different  shapes,  but  they  differ  by  an  order  of  magnitude. 
Unfortunately  at  this  time  we  are  unable  to  verify  the  NAD  terms  from  Brooklyn  by  our 
own  calculations  in  Columbus,  which  was  to  be  the  first  exploration  in  finding  the  error 
in  previous  calculations.  The  negative  result  from  Columbus  is  probably  a  result  of 
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running  it  with  faulty  input  files.  However,  correction  of  the  error  in  the  second 
exploration  (the  conversion  of  Cartesian  to  Jacobi  derivative  coupling  terms)  produced 
results  that  agreed  with  theory;  thus  the  NAD  terms  from  Brooklyn  are  demonstrably 
correct. 

The  Jacobi  Transformation 

Let  us  now  examine  the  NAD  terms  in  Jacobi  coordinates.  With  six  Cartesian 
coordinates  we  are  able  to  define  any  possible  configuration  of  B+Hi  within  the  plane. 
When  creating  inputs  for  Columbus  and  Brooklyn  in  Figure  5  we  were  able  to  restrict  the 
placement  of  the  system  in  the  plane  to  minimize  the  number  of  Cartesian  coordinates 
required  to  specify  an  internal  geometry.  This  restriction  was  acceptable  for  the  purposes 
of  simplifying  inputs  because  the  potential  energy  surfaces  and  NAD  terms  only  depend 
on  internal  geometry;  that  is,  they  only  depend  on  the  relative  positions  of  the  atoms  to 
each  other,  not  where  they  are  with  respect  to  the  origin. 

In  order  to  convert  from  one  coordinate  system  to  another,  we  need  to  lift  the 
restriction  of  how  the  system  is  oriented  with  respect  to  the  origin.  .  The  usual  Jacobi 
coordinates  are  internal;  they  orient  the  atoms  with  respect  to  each  other  but  specify  no 
information  whatsoever  about  where  in  the  plane  the  system  is  found.  We  need  to 
augment  our  three  Jacobi  coordinates  with  an  overall  displacement  and  an  overall  rotation 
We  choose  the  center  of  mass  of  the  system  as  the  point  to  measure  displacement, 
although  any  point  that  is  fixed  in  the  center  of  mass  frame  of  reference  will  yield  the 
same  result.  At  this  point  the  choice  of  rotation  angle  reference,  0,  is  arbitrary  inasmuch 
as  it  is  independent  of  the  other  five  Jacobi  coordinates.  The  choice  of  Dr.  Yarkony  is 
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Figure  15.  The  Six  Jacobi  Coordinates 

illustrated  in  Figure  15.  His  co  is  an  average  of  9j  and  02  (Weeks,  2004).  This 
arrangement  also  allows  y  to  be  alternatively  defined  as  02-01-  Thus  the  six  Jacobi 
coordinates  are:  R,  r,  y,  xcra,  ycm,  and  co.  These  Jacobi  coordinates  are  independent;  that  is, 
any  one  can  be  changed  without  forcing  any  of  the  others  to  change.  For  the  sake  of 
being  arbitrary,  one  could  just  as  easily  have  chosen  either  9i  or  92  as  the  reference  of 
rotation,  as  either  one  is  also  independent  of  the  other  five  coordinates.  But  as  we  will 
show,  in  order  to  separate  the  internal  and  external  Hamiltonians,  the  choice  of  reference 
angle  must  be  specifically  co  =  0i ;  this  was  the  cause  of  error  in  the  coupling  angles 
calculated  from  the  Jacobi  NAD  terms.  The  six  Cartesian  coordinates  are:  xB  and  yB,  the 
position  of  boron;  xHi  and  yHi  the  position  of  one  of  the  hydrogens;  and  xH2  and  yH2  the 
position  of  the  other  hydrogen. 
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The  output  of  Brooklyn  or  Columbus  gives  us  the  Cartesian  NAD  terms;  but 
mixing  them  all  to  come  up  with  Jacobi  derivatives  is  a  formidable  task.  A  Jacobi 
derivative  coupling  term  is  extracted  from  the  Cartesian  derivative  coupling  terms  by: 


dxB  d 
dR  dx„ 


dR  dyB 


dR  dx. 


dR  dyh 


dR  dxu 


dR  dyh 


(42) 


where  R  can  represent  any  of  the  Jacobi  coordinates.  Inspection  of  equation  (42)  reveals 
that  this  equality  arises  from  the  chain  rule  for  partial  derivatives.  Equation  (43)  reveals 
that  this  is  equal  to  a  product  of  Cartesian  derivative  coupling  terms  and  derivatives  of  the 
Cartesian  coordinates  with  respect  to  the  Jacobi  coordinate  in  question.  Given  the  vectors 


(44) 


and 

r  =  (XH2  —  XHl-  y H2  ~  y HI  )’  (45) 

the  Jacobi  coordinates  in  terms  of  the  Cartesians  are 

R  =  \r\,  (46) 

r  =  Irl ,  (47) 
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It  is  very  difficult  if  not  impossible  to  analytically  solve  for  the  derivatives  in  equation 
(43)  because  the  above  transformation  is  not  invertible.  Fortunately  the  reverse  operation, 
getting  the  derivatives  of  the  Jacobi  coordinates  with  respect  to  the  Cartesian  coordinates 
is  fairly  straight  forward.  The  matrix  that  contains  all  the  partial  derivatives  of  one  set  of 
coordinates  with  respect  to  another  set  is  called  the  Jacobian.  The  Jacobian  for  the 
Cartesian-to-Jacobi  operation  looks  like 
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This  matrix  dotted  into  a  vector  containing  the  Cartesian  derivative  coupling  terms  leads 
to  six  values  that  have  the  form  of  the  right-hand-side  of  equation  (43): 


dxB 

dyB 

dr 

dr 

dxB 

dyB 

dR 

dR 

(Weeks,  2004) 

where  the  dR  /  dx:  are  known  analytically  from  the  R=R(x;)  transformation  in  equations 


(46)-(51).  Careful  consideration  shows  that  the  product  of  these  matrices  is  indeed  the 
identity  matrix,  which  confirms  that  going  from  Cartesian  to  Jacobi  to  Cartesian  or  vice 
versa  is  an  identity  operation  and  should  leave  the  coordinates  unchanged.  It  then 
becomes  a  simple  matter  of  linear  algebra  to  use  this  new  Jacobian  in  conjunction  with 
the  Cartesian  NAD  terms  to  come  up  with  the  Jacobi  NAD  terms  as  per  equation  (43). 
Figure  16  shows  the  NAD  terms  from  Figure  13  (those  calculated  by  Columbus)  after 
having  been  converted  to  Jacobi  coordinates.  Notice  how  erratic  the  angular  derivatives 
are  (d/dy  and  d/dro).  Figure  17  shows  the  NAD  terms  from  Figure  14  (those  calculated  by 
Brooklyn )  after  having  been  converted  to  Jacobi  coordinates.  Note  that  all  the  NAD 
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Figure  16.  Jacobi  NAD  Terms  from  Columbus 
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Figure  17.  Jacobi  NAD  Terms  from  Brooklyn 


39 


terms  change  smoothly  across  the  energy  curve  slice,  especially  the  angular  terms.  Note 
also  the  much  smaller  scale  that  the  Brooklyn  NAD  terms  take  on.  Keep  in  mind  that 
these  Jacobi  NAD  terms  were  calculated  with  Dr.  Yarkony’s  co  from  Figure  15,  so  theyare 
not  necessarily  correct  (specifically  the  d/dy  term).  Initial  calculation  using  these  values 
led  to  the  incorrect  coupling  angle  as  the  next  section  will  illustrate. 

Calculating  the  Coupling  Angle 

The  derivative  coupling  terms  we  have  calculated  are  used  to  create  the  coupling 
angle  0(R)  in  equation  (37),  which  is  then  used  to  form  the  rotation  matrix  in  equation 
(24).  This  rotation  matrix  then  diagonalizes  the  kinetic  energy  operator  and  mixes  the 
potential  energy  surfaces  in  the  potential  energy  operator,  placing  them  in  the  diabatic 
representation.  Because  the  integral  in  equation  (37)  is  path  independent,  we  can  choose 
to  integrate  from  an  arbitrary  starting  point  through  the  derivative  coupling  terms  to  the 
point  (R,  r,  y)  in  straight  lines  parallel  to  the  axes  (See  Figure  18).  This  line  integral  can 


d/dR 


Figure  18.  Path  Integral  to  Calculate  0(R)  in  NAD-Term  Space 
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be  made  even  simpler  by  choosing  to  place  the  starting  point  for  the  integration  in  the 
plane  y  =  ti/2.  Notice  in  Figure  17  that  the  NAD  terms  with  respect  to  r  and  R  are  nearly 
zero  at  y  =  nil.  This  is  true  for  all  values  of  R  and  r.  Thus  the  integral  from  0(Ro)  to  A 
and  A  to  B  is  negligible,  and  the  entire  integral  is  equal  to  the  integral  from  B  to  0(R). 
This  implies  that  0(R)  need  only  be  integrated  along  y  (Weeks,  2004). 

In  his  work  on  the  B+FF  system  with  hydrogen  at  its  equilibrium  distance, 
Alexander  calculated  that  whether  the  energy  eigenvalues  of  the  electronic  Hamiltonian 
were  positive  or  negative  would  depend  on  how  close  the  hydrogens  came  to  either  end  of 
the  antisymmetric  boron  atomic  p-orbital  (1993:6018).  Recall  that  these  energy 
eigenvalues  form  the  potential  energy  surfaces  for  the  nuclear  Hamiltonian.  If  the 


y  <  ti/2,  positive  energy  y  >  ti/2,  negative  energy 
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Figure  19.  How  y  Affects  Electronic  Eigenvalues 
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hydrogen  swings  closer  to  the  positive  lobe  there  is  a  repulsion  and  hence  a  positive 
energy.  If  it  swings  closer  to  the  other  lobe  the  opposite  is  true  (see  Figure  19).  Since  the 
potential  energy  surfaces  must  be  continuous,  this  implies  that  in  the  collinear 
configuration  and  in  the  “T”  configuration  all  the  surfaces  must  go  through  zero. 
Alexander  also  derived  that  the  off-diagonal  (coupling)  potential  energy  surface  would  be 
related  to  the  coupling  angle  by: 

El2  =sin{®{R ))cos(Q(R ))-(El  -  £2  )(1993:6020).  (55) 

In  order  for  the  surface  to  go  to  zero,  this  implies  that  0(R)  must  go  to  0  or  7t/2  at  y  =  n/2 
and  y  =  0.  Because  we  put  the  arbitrary  starting  point  of  the  integral  at  y  =  7i/2,  0(R,  r, 
7i/2)  will  always  be  equal  to  zero.  Figure  20  shows  Alexander’s  calculation  of  coupling 
angle  as  a  function  of  Jacobi  angle.  For  below  values  of  about  R  =  6.5  the  line  integral 
goes  back  to  zero  and  for  values  above  it  goes  to  tt/2,  as  shown  in  Figure  20.  The  y  NAD 


y (degree) 

Figure  20.  Alexander's  Calculation  of  the  Coupling  Angle  (1993:6019) 
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term  as  calculated  in  Brooklyn  for  the  R=5,  r=  1 .402  slice  of  the  energy  surface  is  shown 
in  Figure  21  (this  is  extracted  from  Figure  17).  For  most  of  the  slice  the  NAD  term  is 
negative,  with  only  a  small  portion  of  the  curve  near  zero  being  positive.  The  result  of 


<1A'  |~  |2A'> 


Figure  21.  y  NAD  Terms  from  Brooklyn 


e  (ud) 


Figure  22.  Coupling  Angle  from  Brooklyn  NAD  Terms 

the  line  integral  across  this  slice  is  shown  in  Figure  22;  the  value  approaches  neither  0  nor 
n/2  (this  plot  is  actually  the  absolute  value  of  the  line  integral  for  the  sake  of  comparing 
with  Alexander’s  data).  This  disagreement  with  theory  implies  that  either  the  derivative 
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coupling  terms  produced  by  Brooklyn  are  wrong  or  the  method  of  transforming  them  to 
Jacobi  coordinates  is  wrong. 

Origin  Dependency 

Recall  that  when  we  oriented  the  system  in  the  plane  with  Jacobi  coordinates,  we 
specified  the  angle  co  =  (0i+02)/2,  but  there  were  several  choices  of  angles  that  were 
independent  of  the  other  five  Jacobi  coordinates.  Figure  23  shows  the  Jacobi  NAD  terms 
calculated  with  co  =  0 1 .  Compare  this  figure  with  Figure  17,  which  was  calculated  with  co 
=  (0i+02)/2.  All  the  derivative  coupling  terms  are  exactly  the  same  with  the  exception  of 
the  y  term.  Even  so  the  y  term  appears  to  have  the  same  line  shape,  but  the  terms  are 
offset  by  a  common  constant  (0.5).  Clearly  by  choosing  different  but  equally  valid 


-------  d/dr 
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Figure  23.  Jacobi  NAD  Terms  from  Brooklyn,  co  =  0j 
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coordinate  systems  we  affect  some  of  the  derivative  coupling  terms.  This  contingency  on 
the  choice  of  coordinates  is  known  as  an  origin  dependency. 

Origin  dependency  is  best  understood  by  the  following  example.  Consider  a 
linear  system  of  two  masses  as  in  Figure  24.  We  can  construct  two  separate  coordinate 


systems  to  describe  their  positions.  The  Cartesian  coordinate  set  describes  the  system  as 
(xi,  X2)  where  xj  is  the  distance  of  mass  mi  from  the  origin  and  X2  is  the  distance  of  mass 
m2  from  the  origin.  The  second  system,  which  might  be  dubbed  the  center-of-mass 
coordinate  set  describes  the  system  as  (r,  s)  where  r  is  the  distance  between  the  masses 
and  s  is  the  distance  from  the  origin  to  some  point  relative  to  the  positions  of  the  two 
masses.  If  s  is  chosen  such  that  s  =  mi*xi+m2*X2  then  s  is  the  distance  to  the  center  of 
mass,  but  in  general  the  center-of-mass  set  is  related  to  the  Cartesian  set  by: 

r  =  x2-xl  (56) 

and 
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s  =  a*xl+b*x2 

If  we  express  this  relationship  in  matrix  form, 

fa  bf  f  x, } 


VrJ  f-l  1 )  \x2) 

we  can  easily  invert  it  in  order  to  find  the  Cartesian  set  in  terms  of  the  center-of-mass  set: 

(xA  1  ( 1  -b']  M 


From  this  we  can  extract  the  Jacobian: 


x2J  a  +  6^1  a  J 


d.v'[ 

dx,A 

f  1 

1 

ds 

ds  _ 

a  +  b 

a  +  b 

dx1 

dx2 

-b 

a 

dr 

dr  J 

Ka+b 

a+b 

Jacobian  = 


Now  suppose  we  have  the  derivative  of  some  function,  vp  with  respect  to  the  Cartesian 
coordinates,  but  we  want  its  derivative  with  respect  to  the  center-of-mass  coordinates. 
We  just  dot  the  Jacobian  into  the  derivatives  (compare  equation  (43)): 

fd^ 

ds  _  ds  ds  #  ^xi  (, 

dT*  d.Vj  c)x2  dvF 

V dr  )  V  dr  dr  J  ^dx2y 

The  derivatives  evaluate  to 

dW  1  dT*  1  d¥ 

- = - + -  (i 

d s  a  +  b  dx,  a  +  b  dx2 


dW  _  -b  d¥  |  a  d¥ 
dr  a  +  b  dxx  a  +  b  dx2 
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As  long  as  the  sum  of  a  and  b  are  the  same,  the  derivative  with  respect  to  s,  the  external 
coordinate,  remains  the  same;  but  the  derivative  with  respect  to  r,  the  internal  coordinate, 
varies  as  the  ratio  of  a  to  b  varies. 

This  simple  example  helps  to  explain  the  relationship  between  the  co  and  the  y 
NAD  terms.  The  internal  angular  coordinate,  y,  was  defined  as  02-0 1  whereas  the  external 
angular  coordinate,  co,  was  defined  as  a*0i+b*02.  As  long  as  a+b  =  1  the  co  NAD  terms 
remained  the  same,  while  the  y  NAD  terms  varied  with  the  ratio  of  a  to  b.  This  implies 
that  the  choice  of  a  and  b  is  not  arbitrary  when  calculating  the  y  NAD  term.  Part  of  the 
reason  that  co  =  (0i+02)/2  was  not  discounted  immediately  upon  finding  the  error  in  the 
coupling  angle  was  that  the  value  of  the  d/dco  term,  which  is  a  measure  of  angular 
momentum,  was  compared  to  and  found  matching  the  angular  momentum  calculated 
using  other  methods  and  it  was  assumed  that  all  other  NAD  terms  must  also  be  correct. 

The  correct  choice  of  a  and  b  to  extract  the  correct  derivative  coupling  terms  and 
consequently  the  correct  coupling  angle  is  buried  in  another  form  of  the  Hamiltonian 
which  is  generalized  for  use  in  three  dimensions: 
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(64) 


where  the  first  term  is  the  linear  momentum  of  the  boron  with  respect  to  the  hydrogen,  the 
second  term  is  the  vibrational  motion  of  the  hydrogen  molecule,  the  third  is  the  angular 
momentum  of  the  hydrogen,  the  fourth  is  the  angular  momentum  of  the  system,  the  fifth 
is  the  potential  due  to  the  electrons,  and  the  sixth  is  a  spin-orbit  term  (Niday,  1999:20). 
The  third  term  in  the  Hamiltonian  is  associated  with  y  and  the  fourth  term  is  associated 
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with  co.  Coming  to  this  form  of  the  Hamiltonian  in  which  the  internal  and  external 
angular  momenta  are  separated  required  using  the  Euler  angles  as  part  of  the  coordinate 
system.  Since  we  are  restricting  our  system  to  a  plane  we  need  not  use  all  of  the  Euler 
angles;  however,  we  still  need  one  of  them,  p,  to  account  for  the  system’s  rotation  in  the 
plane.  The  Euler  angle  P  is  the  angle  0i  from  Figure  15.  Thus  in  order  to  successfully 
separate  the  Hamiltonian  to  arrive  at  the  correct  coupling  angles,  we  must  set  a  to  1  and  b 
to  0;  that  is,  co  =  0j.  Figure  25  shows  the  y  NAD  terms  for  R  =  5,  r  =  1.402,  and  the 
resultant  coupling  angle  after  setting  co  to  equal  0i.  In  agreement  with  Millard  Alexander’s 


Figure  25.  y  NAD  Terms  and  Line  Integral  After  Effecting  Rotation 

data,  the  coupling  angle  now  goes  to  zero  at  y  =  0.  Figure  26  shows  the  line  integrals  after 
rotation  for  a  variety  of  lengths  of  R.  Table  1  compares  their  value  at  y  =  0  to  0  or  7i/2. 
Since  these  coupling  angles  behave  almost  exactly  as  Alexander  calculated,  we  can 
conclude  that  the  derivative  coupling  terms  given  by  Brooklyn  are  correct. 
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Table  1.  Line  Integrals  Expected  vs.  Actual  Values 


R 

Expected  (rad) 

Actual  (rad) 

Difference  (rad) 

3 

0 

0.028819 

-0.02882 

4 

0 

0.006078 

-0.00608 

5 

0 

0.012663 

-0.01266 

6 

0 

0.04257 

-0.04257 

7 

1.570796327 

1.588636 

-0.01784 

8 

1.570796327 

1.619886 

-0.04909 

9 

1.570796327 

1 .595284 

-0.02449 

10 

1.570796327 

1.587158 

-0.01636 

Calculating  the  Non-Adiabatic  Surfaces 

Having  verified  the  coupling  angles  we  may  now  proceed  to  construct  the  non- 


adiabatic  surfaces  which  we  have  sought.  The  coupling  angle,  0(R)  is  used  to  construct 
the  rotation  matrix  in  equation  (24).  The  rotation  matrix  is  then  applied  to  both  the 
kinetic  and  potential  energy  operators  as  in  equation  (25).  Rotation  of  the  potential 
energy  operator  creates  the  diabatic  potential  energy  operator: 
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.(65) 


The  new  diagonal  terms  are  the  diabatic  energy  surfaces,  and  only  the  1A”  surface  which 
did  not  mix  remains  the  same.  The  off-diagonal  terms  are  the  coupling  surface.  Although 
not  a  proper  surface,  the  coupling  surface  gives  a  measure  of  the  likelihood  that  the  1A’ 
and  2A’  surfaces  will  couple;  that  is,  that  a  wave  function  found  on  one  will  work  its  way 
onto  the  other. 
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Application  of  this  method  to  the  slices  shown  in  Figure  1 1  results  in  the  surfaces 


shown  in  Figure  27.  This  figure  includes  only  the  1A’  and  2A’  surfaces  since  the  1A” 
surface  does  not  mix  and  is  left  as  before  in  its  adiabatic  state.  Once  again  these  are  two- 


Figure  27.  Surface  and  Contour  Plots  of  the  1A'  and  2 A'  Diabatic  Surfaces  and  the 


Coupling  Surface 
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dimensional  slices  of  three-dimensional  surfaces,  the  slice  being  taken  at  r  =  1.402  bohr. 
The  first  column  has  the  surface  plots  and  the  second  column  has  overhead  contour  plots. 
The  1  A’  and  2 A’  surfaces  are  plotted  on  the  range  of  R  =  5  to  10  and  y  =  0  to  ji.  To  this 
figure  we  have  also  added  the  feature  we  had  been  missing  from  the  adiabatic  picture:  the 
coupling  surface.  The  coupling  surface  shows  more  coupling  as  R  gets  smaller  and  less 
as  R  gets  larger.  It  also  shows  that  at  y  =  ti/2  there  is  no  coupling.  This  is  significant 
because  this  is  the  region  on  the  1A’  surface  where  the  potential  well  is  deepest  and  the 
Van  der  Waals  molecule,  should  it  be  formed,  is  most  likely  to  be  found.  That  implies 
that  there  is  very  little  possibility  of  the  molecule  “leaking”  out  onto  the  antibonding  2A’ 
surface,  which  would  seem  to  confirm  that  B+FT  is  a  good  candidate  for  being  a  HEDM. 
However,  should  the  hydrogen  molecule  rotate  to  the  left  or  right  with  respect  to  the 
boron,  the  coupling  becomes  very  large. 

These  surfaces  can  be  compared  to  the  surfaces  that  Alexander  calculated  and  be 
found  matching  (1993;  6021-24).  This  too  is  significant  because  it  means  that  the  NAD 
terms  calculated  in  Brooklyn  by  Dr.  Yarkony  are  correct.  Furthermore  this  verifies  that 
the  definition  of  to  as  the  Euler  angle  P  is  correct  and  produces  the  correct  coupling  angle. 
This  means  that  the  balance  of  the  surfaces  for  lengths  of  r  different  from  the  equilibrium 
distance  can  be  mapped  out  with  confidence  using  this  method. 
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IV.  Conclusion 


Adiabatic  chemistry  has  been  the  staple  of  quantum  chemistry  for  several  decades. 
The  Born-Oppenheimer  approximation  has  made  problems  tractable  on  computers. 
Nevertheless,  as  computers  become  more  capable  and  we  demand  more  accuracy  in  our 
calculations  of  molecular  wave  functions  and  energy  surfaces,  diabatic  chemistry  or  at 
least  adiabatic/diabatic  hybrid  chemistry  is  becoming  increasingly  important.  As  it  works 
its  way  into  code  and  models  we  must  be  aware  of  the  various  problems  associated  with 
its  implementation,  like  this  origin  dependency.  Extreme  care  must  be  taken  when 
calculating  the  coupling  between  non-adiabatic  energy  surfaces  so  that  this  origin 
dependency  or  one  similar  to  it  does  not  creep  its  way  into  the  results  and  papers  of 
unaware  researchers. 

As  for  this  project,  we  can  conclude  that  the  derivative  coupling  terms  created  by 
Brooklyn  are  indeed  correct.  We  have  calculated  a  subset  of  the  coupling  angles  and 
diabatic  surfaces  from  these  data  and  found  them  agreeing  with  the  work  of  Millard 
Alexander.  Follow-on  projects  will  now  be  able  to  use  the  data  to  complete  the  full  three- 
dimensional  surfaces.  These  surfaces  can  then  be  used  to  construct  the  scattering  matrix 
and  characterize  whether  B+FT  can  be  used  as  a  high  energy  density  fuel. 

On  the  subject  of  Columbus  it  is  unfortunate  that  we  were  not  able  to  create  our 
own  derivative  coupling  terms  similar  to  those  created  by  Brooklyn.  Many  man-hours 
were  used  to  tweak  Columbus’  switches  and  knobs  and,  although  progress  was  made,  the 
right  inputs  could  not  be  found.  Hopefully  with  more  contact  and  communication  with 
the  maintainers  of  Columbus  future  researchers  in  the  Department  of  Physics  in  need  of 
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derivative  coupling  terms  will  be  able  to  master  the  code  and  create  them  according  to 
their  needs. 
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Appendix  A:  The  Scattering  Matrix 


The  end  of  finding  the  non-adiabatic  potential  surfaces  is  making  the  scattering 
matrix.  Because  the  scattering  operator  exists  in  linear  momentum  space  crossed  with 
linear  momentum  space,  its  matrix  representation  is  dense  and  infinitely  large.  We  can 
therefore  only  calculate  a  certain  portion  of  the  matrix  around  values  we  find  most 
interesting.  Further,  we  can  only  calculate  the  matrix  in  discrete  steps  and  must 
interpolate  the  other  values  if  we  should  need  them. 

The  scattering  problem  begins  by  specifying  |  )  and  |  'Poltf ) ,  the  wave  packets 

of  the  initial  and  final  states,  respectively  (Niday,  1999:6).  These  states  must  then  be 
propagated  to  their  asymptotic  limits,  reactants  to  t  =  -°°  and  products  to  t  =  °°  using  the 
asymptotic  Hamiltonian: 

|  ^Intermedia, e  \>  =  £xp  {_  !^}|  XJ/  }  (66) 

and 

|f  "“)  =  (67) 

These  intermediate  states  are  then  taken  from  infinity  and  propagated  toward  interaction 
via  the  Hamiltonian  created  with  the  potential  surfaces: 

TM  =  exp  {- — }|  \  (68) 

fi  1  ' 

and 

I'P  >  =  exp {- — }|  \  .  (69) 

fi  1  ' 
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The  scattering  matrix  element  is  a  “bridge”  between  these  two  states,  so  in  order  to 
calculate  it,  we  must  calculate  a  correlation  between  the  two  states: 


Cfr{t)  =  ('¥l 


exP{-“~}  ‘ 
n 


(70) 


where  the  reactants  and  products  have  now  been  indexed  by  y  and  y  to  differentiate 
specific  beginning  and  ending  states.  The  scattering  matrix  operator  element  is  then 
calculated  as: 
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(71) 


where  E  is  the  total  energy,  ky  is  the  relative  momentum  between  reactants,  ky  is  the 
relative  momentum  between  products,  7+  (±  ky )  are  expansion  coefficients  used  to 

construct  the  initial  states,  and  py  are  the  reduced  masses. 
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Appendix  B:  The  Many-Body  Problem 

Although  we  have  traced  the  theory  of  non-adiabatic  energy  surfaces  from 
beginning  to  end,  we  are  unable  to  put  it  to  use  unless  we  actually  calculate  the  electronic 
wave  functions  themselves. 

The  student  of  quantum  mechanics  is  well  acquainted  with  the  derivation  of 
hydrogenic  orbitals.  He  has  also  been  reminded  time  and  time  again  by  his  professors 
that  this  two-body  problem  is  the  only  one  of  its  kind  for  which  one  can  solve  analytically 
the  molecular  wave  functions.  Even  the  addition  of  a  single  electron  puts  calculation  of 
molecular  wave  functions  squarely  in  the  hands  of  the  computer.  In  the  B+H2  system  we 
have  10  bodies  including  nuclei  and  electrons.  Fortuitously  we  are  not  without  help — a 
number  of  methods  to  iteratively  approximate  molecular  wave  functions  are  available. 
Hartree-Foek 

The  most  basic  way  to  determine  molecular  wave  functions,  and  the  method  upon 
which  more  complicated  models  are  based,  is  the  Hartree-Foek  method  (Szabo, 
1996:108-230).  This  method  solves  the  Hartree-Foek  eigenvalue  equation, 

-ivf  -X—  +v"r(i)\x(x,)  =  ez{x,)  (72) 

v  ^  a  ria  j 

Tip  j-U 

where  the  v  is  the  average  potential  experienced  by  the  i  electron  due  to  the  other 
electrons.  The  quantity  in  parentheses  is  known  as  the  Fock  operator.  By  introducing 
this  average  interaction,  the  energy  s  is  merely  a  sum  of  the  individual  orbital  energies  for 
each  electron,  and  the  wave  function,  XF,  is  a  product  of  the  hydrogenic  type  orbital  of 
each  electron: 
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,  X2-Xn  )  =  Zl  (Xl  )Z2  (X2  )-Zn  (*„  ) 


(73) 


where  x>  represents  a  spin  orbital  and  x,  are  the  coordinates  of  a  given  electron.  A  spin 
orbital  is  the  product  of  an  atomic  orbital  and  a  spin  function,  which  can  be  either  up 
(a(co))  or  down  (P(a>)).  This  form  of  the  wave  function  is  not  yet  complete;  it  does  not 
necessarily  conform  to  the  antisymmetry  principle.  This  principle,  which  is  born  out  of 
the  Pauli  Exclusion  Principle,  states  that 

...v  ...I-  (74) 


In  order  that  obeys  this  principle,  we  introduce  Slater  Determinants,  so  that 


Xit(xt,x2...xn ) 


Zi(xi) 

Zi{x2) 


Z2(x  i) 
Z 2  (X2  ) 


(75) 


These  Slater  Determinants  are  often  written  in  simplified  form: 

v¥(xl,x2...xn)=jz,Z2-Zn)  (7  6) 

where  the  normalization  is  implied.  Hartree-Fock  uses  a  single  Slater  Determinant  to 
determine  the  ground  state  energy.  The  variation  principle  states  that  any  energy  we  find 
by  adding  and  subtracting  atomic  orbitals  will  be  an  upper  bound  for  the  true  energy  of 
the  configuration;  that  is,  the  best  approximation  to  the  molecular  wave  function  is  the 
one  with  the  lowest  energy.  This  can  be  done  iteratively  with  the  help  of  a  software 
package.  An  initial  guess  at  the  linear  combination  of  atomic  spin  orbitals  that  gives  the 
Slater  Determinant  will  give  the  average  field  term  in  the  Fock  operator.  This  Fock 
operator  will  have  eigenfunctions  which  are  Slater  Determinants  different  from  those 
initially  guessed.  These  new  combinations  can  be  used  to  find  a  new  average-field  term 
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again,  and  so  on  until  the  eigenfunctions  are  consistent  with  the  Fock  operator.  Because 
the  object  is  to  iteratively  vary  the  Slater  Determinants  until  the  electric  field  calculated 
from  them  is  consistent  with  the  Fock  operator  whose  eigenfunctions  they  are,  this 
process  is  called  self-consistent  field,  or  SCF.  This  method  is  an  approximation — but  the 
more  atomic  orbitals  that  are  used  the  lower  the  energy  will  be  and  the  better  the 
approximation. 

Configuration  Interaction  (Cl)  and  MCSCF 

The  Configuration  Interaction  (Cl)  method  is  a  refinement  beyond  Hartree-Fock 
that  is  considered  a  deeper  level  of  theory  by  using  multiple  Slater  Determinants  rather 
than  just  one  as  in  Hartree-Fock  (Szabo,  1996:231-270).  Given  a  complete  set  of  spin 
orbitals,  xi(xi),  one  can  construct  any  arbitrary  function  within  the  space: 

*(*,)=!«, *,&)  (77) 

i 

It  can  be  shown  as  well  for  two  variables: 

cIiU'i ,  *2 )  =  X  ai.jXi  (*i  )Zj  (x2 )  (78) 

i.j 

As  before  if  this  is  a  wave  function  we  would  like  it  to  be  antisymmetric,  so  we  use  a 
Slater  Determinant  instead: 

<S>{X],x2)  =  Yjauj\xiX])  (79) 

i<j 


And  in  general 

<S>(xltX2...ZH)=  Zau-nlzj 'j-Xn)  (80) 

i<j...n 

Thus  Cl  goes  beyond  Hartree-Fock  by  defining  the  wave  function  not  only  by  the  single 
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Slater  Determinant  giving  the  lowest  energy,  but  as  a  linear  combination  of  all  possible 
Slater  Determinants  of  n  electrons.  This  yields  a  lower  and  thus  more  accurate  energy  for 
the  configuration  by  the  variation  principle.  If  we  define  |d>0)  as  the  Hartree-Fock 
ground  state  Slater  Determinant,  we  can  express  other  Slater  Determinants  as  excitations 
out  of  this  one.  For  instance,  if  |<I>0)  =  |  Xi—Xi—Xn)  then  we  can 

create  |<b/  J  =  | X\—Xj—Xn ) »  a  singly  excited  determinant,  by  replacing  the  spin  orbital  y, 

with  %j.  We  can  also  create  and  name  doubly,  triply,  quadruply,  et  c.  excited  determinants 
and  define  the  wave  function  as: 

ar  a<b,r<s 

Once  again  we  are  limited  by  the  fact  that  we  do  not  have  an  infinite  number  of  atomic 
orbitals  with  which  to  construct  molecular  wave  functions,  nor  can  we  sum  several 
multiply-excited  determinants  before  we  get  bogged  down.  As  in  the  Hartree-Fock 
method,  a  truncated  basis  of  atomic  orbitals  must  be  chosen,  and  a  limit  on  excited 
determinants  must  be  chosen  as  well.  Many  software  packages  limit  the  sum  to  doubly 
excited  determinants.  Then  the  determinants  as  well  as  the  coefficients  are  varied  until  a 
lowest  upper  bound  is  found  for  the  energy.  This  process  of  optimization  is  called 
multiconfiguration  self-consistent  field,  or  MCSCF. 
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Appendix  C:  Columbus’  Calculation  of  the  Derivative  Coupling  Term 

Lischka  et  al.  describe  the  method  Columbus  uses  to  calculate  the  NAD  terms 
which,  while  more  in  depth,  takes  a  similar  form  at  least  in  part  to  the  method  used  by 
Koppel  (2004:7323).  From  Cl  we  know  that  the  electronic  wave  functions  are  linear 
combinations  of  other  functions.  In  Appendix  B  they  were  referred  to  as  various  Slater 
Determinants.  In  the  text  of  Lischka  et  al.  the  term  configuration  state  function  (CSF)  is 
used  (2004:7322).  Thus  the  electronic  wave  functions  take  the  form 

k)=Xc«K)  (82) 

n 

where  vpn  is  the  CSF.  Substituting  into  the  form  of  the  NAD  term  we  get: 

p„  =  (83) 

n  \  m  J 

Thus  the  NAD  term  can  be  further  broken  down: 

Pij  =  E  c»  E  v(c'»  ) + E  c»  E cl™  I  ^  v™ )  •  (84) 

n  m  n  m 

Since  the  vp’ s  are  orthonormal,  the  first  term  is  just  an  inner  product  leaving  only  the 
coefficients  so  that 

Z<<^Eri^)|r,„>={c'|vcJ).  (85) 

n  m 

Lischka  et  al.  use  the  same  rationale  as  Koppel  to  rewrite  this  term  as 

■  (86) 

Thus  the  NAD  term  can  be  written  as  two  distinct  terms: 
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(87) 


ptj  = 


—  l 

c 


VH  c- 


EJ~E > 


+Zc-c™(^»  lv^) 


the  former  being  called  the  Cl  term  and  the  latter  the  CSF  term  (Lischka  2004:7323). 

This  is  how  the  NAD  term  is  calculated  in  Columbus  and  most  likely  in  Brooklyn  as  well. 
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